# devtools::install_github("cran/drc", force = TRUE)
# devtools::install_github("cran/alr3", force = TRUE)
# devtools::install_github("cran/NRAIA", force = TRUE)
# devtools::install_github("cran/nlrwr", force = TRUE)
library(nlrwr)
Loading required package: alr3
Loading required package: car
Loading required package: carData
Loading required package: drc
Loading required package: MASS
Attaching package: ‘MASS’
The following object is masked from ‘package:alr3’:
forbes
'drc' has been loaded.
Please cite R and 'drc' if used for a publication,
for references type 'citation()' and 'citation('drc')'.
Attaching package: ‘drc’
The following objects are masked from ‘package:stats’:
gaussian, getInitial
Loading required package: HydroMe
Loading required package: lattice
Loading required package: lmtest
Loading required package: zoo
Attaching package: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
Loading required package: NISTnls
Loading required package: nlme
Loading required package: nls2
Loading required package: proto
Loading required package: nlstools
'nlstools' has been loaded.
IMPORTANT NOTICE: Most nonlinear regression models and data set examples
related to predictive microbiolgy have been moved to the package 'nlsMicrobio'
Loading required package: NRAIA
Registered S3 method overwritten by 'NRAIA':
method from
plot.profile.nls stats
Attaching package: ‘NRAIA’
The following object is masked from ‘package:nlstools’:
plotfit
Loading required package: sandwich
'nlrwr' has been loaded.
Please cite R and 'nlrwr' if used for a publication,
for references type 'citation()' and 'citation('nlrwr')'.
Attaching package: ‘nlrwr’
The following object is masked from ‘package:nlstools’:
confint2
## Fig. 1.1.
plot(num.fish ~ spawn.biomass,
data = M.merluccius,
xlab = "Spawning biomass (1000 tonnes)",
ylab = "Recruitment (million fish)")
## Fig. 1.2.
plot(biomass ~ x,
data = RScompetition,
log = "",
xlab = Density ~ (plants/m^2),
ylab = Biomass ~ of ~ sensitive ~ biotype ~ (g/plant),
pch = as.numeric(as.factor(RScompetition$z)))
## Fig. 1.3.
xyplot(DryMatter ~ Dose | Herbicide,
data = S.alba,
scales = list(x = list(log = TRUE)),
ylab = "Dry matter (g/pot)",
xlab = "Dose (g/ha)")
L.minor
## Fig. 2.1.
plot(rate ~ conc, data = L.minor,
ylab = "Uptake rate (weight/h)",
xlab = Substrate ~ concentration ~ (mmol ~ m^-3))
L.minor.m1 <- nls(rate ~ Vm * conc/(K + conc),
data = L.minor,
start = list(K = 20, Vm = 120),
trace = TRUE)
624.3282 (1.32e+00): par = (20 120)
244.5460 (2.24e-01): par = (15.92382 124.5715)
234.5198 (2.91e-02): par = (17.25299 126.4388)
234.3595 (5.72e-03): par = (17.04442 125.9618)
234.3533 (1.11e-03): par = (17.08574 126.0467)
234.3531 (2.15e-04): par = (17.07774 126.0302)
234.3531 (4.19e-05): par = (17.0793 126.0334)
234.3531 (8.14e-06): par = (17.07899 126.0328)
deviance(L.minor.m1)
[1] 234.3531
d2 <- sum((L.minor$rate -predict(L.minor.m1))^2)
all.equal(deviance(L.minor.m1), d2)
[1] TRUE
logLik(L.minor.m1)
'log Lik.' -24.86106 (df=3)
coef(L.minor.m1)
K Vm
17.07899 126.03276
summary(L.minor.m1)
Formula: rate ~ Vm * conc/(K + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 17.079 2.953 5.784 0.00117 **
Vm 126.033 7.173 17.570 2.18e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.25 on 6 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 8.141e-06
fitted(L.minor.m1)
[1] 18.06066 28.56474 38.52679 71.09461 78.03806 87.78424 91.62683 116.28685
attr(,"label")
[1] "Fitted values"
concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 10))
predict(L.minor.m1, newdata = data.frame(conc = concVal))
[1] 18.06066 75.09908 92.70509 101.26607 106.32776 109.67157 112.04518
[8] 113.81733 115.19094 116.28685
## Fig. 2.2.
plot(rate ~ conc,
data = L.minor,
ylim = c(10, 130),
ylab = "Uptake rate (weight/h)",
xlab = Substrate ~ concentration ~ (mmol ~ m^-3))
lines(L.minor$conc, fitted(L.minor.m1))
## Fig. 2.3.
plot(rate ~ conc,
data = L.minor,
ylim = c(10, 130),
ylab = "Uptake rate (weight/h)",
xlab = Substrate ~ concentration ~ (mmol ~ m^-3))
concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 100))
lines(concVal, predict(L.minor.m1, newdata = data.frame(conc = concVal)))
abline(h = coef(L.minor.m1)[2], lty = 2)
L.minor.m1con <- nlsContourRSS(L.minor.m1)
0%100%
RSS contour surface array returned
## Fig. 2.4.
par(pty = "s")
plot(L.minor.m1con, col = FALSE, nlev = 10)
## this code should be runned in the clean R session
## library(nlrwr) cause error
# data(L.minor, package = "nlrwr")
# L.minor.m4 <- glm(rate ~ I(1/conc),
# data = L.minor,
# family = gaussian("inverse"))
# summary(L.minor.m4)
## 2.1
y0_start = mean(RGRcurve[RGRcurve$Day == 0, "RGR"])
y1 = mean(RGRcurve[RGRcurve$Day == 1, "RGR"])
b_start = 1 / log(y1/y0_start)
e2.1 <- nls(RGR ~ y0 * exp(Day/b),
data = RGRcurve,
start = list(y0 = y0_start, b = b_start))
plot(RGRcurve$Day, RGRcurve$RGR, xlab = "Day", ylab = "RGR")
nd <- with(RGRcurve, seq(min(Day), max(Day), length.out = 101))
lines(nd, predict(e2.1, newdata = data.frame(Day = nd)), col = 2)
## 2.3
plotfit(L.minor.m1)
## Fig. 3.1.
# install.packages("NISTnls")
data(Chwirut2, package = "NISTnls")
plot(y ~ x,
data = Chwirut2,
xlab = "Metal distance",
ylab = "Ultrasonic response")
## Fig. 3.2.
expFct <- function(x, beta1, beta2, beta3) {
exp(-beta1 * x)/(beta2 + beta3 * x)
}
plot(y ~ x,
data = Chwirut2,
xlab = "Metal distance",
ylab = "Ultrasonic response",
ylim = c(0, 100)
)
curve(expFct(x, beta1 = 1, beta2 = 0.01, beta3 = 1), add = TRUE, col = 2)
## Fig. 3.3.
plot(y ~ x, data = Chwirut2, xlab = "Metal distance",
ylab = "Ultrasonic response", ylim = c(0, 100))
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 1), add = TRUE, lty = 2,
col = 2)
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 0.1), add = TRUE, lty = 3,
col = 3)
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 0.01), add = TRUE, lty = 4,
col = 4)
curve(expFct(x, beta1 = 0.2, beta2 = 0.01, beta3 = 0.01), add = TRUE, lty = 1,
col = 6)
Chwirut2.m1 <- nls(
y ~ expFct(x, beta1, beta2, beta3),
data = Chwirut2,
start = list(beta1 = 0.2, beta2 = 0.01, beta3 = 0.01)
)
summary(Chwirut2.m1)
Formula: y ~ expFct(x, beta1, beta2, beta3)
Parameters:
Estimate Std. Error t value Pr(>|t|)
beta1 0.1665753 0.0383032 4.349 6.56e-05 ***
beta2 0.0051653 0.0006662 7.753 3.54e-10 ***
beta3 0.0121501 0.0015304 7.939 1.81e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.172 on 51 degrees of freedom
Number of iterations to convergence: 5
Achieved convergence tolerance: 4.684e-06
## bad starting example
# Chwirut2.m0 <- nls(
# y ~ expFct(x, beta1, beta2, beta3),
# data = Chwirut2,
# start = list(beta1 = 1, beta2 = 0.01, beta3 = 1))
# summary(Chwirut2.m0)
## nlsLM may solve the bad starting
Chwirut2.m0 <- minpack.lm::nlsLM(
y ~ expFct(x, beta1, beta2, beta3),
data = Chwirut2,
start = list(beta1 = 1, beta2 = 0.01, beta3 = 1))
summary(Chwirut2.m0)
Formula: y ~ expFct(x, beta1, beta2, beta3)
Parameters:
Estimate Std. Error t value Pr(>|t|)
beta1 0.1665765 0.0383033 4.349 6.56e-05 ***
beta2 0.0051653 0.0006662 7.753 3.54e-10 ***
beta3 0.0121500 0.0015304 7.939 1.81e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.172 on 51 degrees of freedom
Number of iterations to convergence: 22
Achieved convergence tolerance: 1.49e-08
grid.Chwirut2 <- expand.grid(beta1 = seq(0.1, 1, by = 0.1),
beta2 = 0.01,
beta3 = seq(0.1, 1, by = 0.1))
grid.Chwirut2
Chwirut2.m2a <- nls2(
y ~ expFct(x, beta1, beta2, beta3),
data = Chwirut2,
start = grid.Chwirut2,
algorithm = "brute-force"
)
Chwirut2.m2a
Nonlinear regression model
model: y ~ expFct(x, beta1, beta2, beta3)
data: Chwirut2
beta1 beta2 beta3
0.10 0.01 0.10
residual sum-of-squares: 60696
Number of iterations to convergence: 100
Achieved convergence tolerance: NA
nls(
y ~ expFct(x, beta1, beta2, beta3),
data = Chwirut2,
start = as.list(coef(Chwirut2.m2a))) # starting value from Chwirut2.m2a
Nonlinear regression model
model: y ~ expFct(x, beta1, beta2, beta3)
data: Chwirut2
beta1 beta2 beta3
0.166578 0.005165 0.012150
residual sum-of-squares: 513
Number of iterations to convergence: 6
Achieved convergence tolerance: 4.429e-06
L.minor.m2 <- nls(rate ~ SSmicmen(conc, Vm, K), data = L.minor)
summary(L.minor.m2)
Formula: rate ~ SSmicmen(conc, Vm, K)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Vm 126.033 7.173 17.570 2.18e-06 ***
K 17.079 2.953 5.784 0.00117 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.25 on 6 degrees of freedom
Number of iterations to convergence: 0
Achieved convergence tolerance: 2.478e-06
expModel <- function(predictor, b, y0) {
y0 * exp(predictor/b)
}
expModelInit <- function(mCall, LHS, data) {
xy <- sortedXyData(mCall[["predictor"]], LHS, data)
lmFit <- lm(log(xy[, "y"]) ~ xy[, "x"])
coefs <- coef(lmFit)
y0 <- exp(coefs[1])
b <- 1/coefs[2]
value <- c(b, y0)
names(value) <- mCall[c("b", "y0")]
value
}
SSexp2 <- selfStart(expModel, expModelInit, c("b", "y0"))
with(RGRcurve, SSexp2(Day, 4, 0.2))
[1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2568051
[8] 0.2568051 0.2568051 0.2568051 0.2568051 0.3297443 0.3297443 0.3297443
[15] 0.3297443 0.3297443 0.3297443 0.4234000 0.4234000 0.4234000 0.4234000
[22] 0.4234000 0.4234000 0.5436564 0.5436564 0.5436564 0.5436564 0.5436564
[29] 0.5436564 0.8963378 0.8963378 0.8963378 0.8963378 0.8963378 0.8963378
[36] 1.4778112 1.4778112 1.4778112 1.4778112 1.4778112 1.4778112
Warning: selfStart caused an ERROR
## ERROR
# getInitial(RGR ~ SSexp2(Day, b, y0), data = RGRcurve)
#RGRcurve.m1 <- nls(RGR ~ SSexp2(Day, b, y0), data = RGRcurve)
#coef(RGRcurve.m1)
More examples on defining self-starter functions for nls() are found in Watkins and Venables (2006) or Venables and Ripley (2002a, pp. 216–217).
MMfct1 <- function(conc, K, Vm) {
numer <- Vm * conc
denom <- K + conc
mean <- numer/denom
partialK <- -numer/(denom^2)
partialVm <- mean/Vm
attr(mean, "gradient") <- cbind(partialK, partialVm)
return(mean)
}
MMfct1(1, 2, 3)
[1] 1
attr(,"gradient")
partialK partialVm
[1,] -0.3333333 0.3333333
L.minor.mgr1 <- nls(rate ~ MMfct1(conc, K, Vm),
data = L.minor,
start = list(K = 1, Vm = 1),
trace = TRUE)
42634.85 (7.03e+00): par = (1 1)
31698.57 (9.47e+00): par = (389.6406 100.9883)
21357.28 (1.23e+01): par = (31.6085 46.75042)
18634.08 (1.03e+01): par = (27.26998 51.07311)
13989.24 (8.02e+00): par = (22.1864 59.6233)
7637.536 (5.64e+00): par = (18.35113 75.51722)
2003.136 (2.75e+00): par = (17.01927 100.4895)
234.3615 (6.31e-03): par = (17.1086 126.0567)
234.3533 (9.56e-04): par = (17.07326 126.021)
234.3531 (1.86e-04): par = (17.08017 126.0352)
234.3531 (3.61e-05): par = (17.07882 126.0324)
234.3531 (7.02e-06): par = (17.07909 126.0329)
summary(L.minor.mgr1)
Formula: rate ~ MMfct1(conc, K, Vm)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 17.079 2.953 5.784 0.00117 **
Vm 126.033 7.173 17.570 2.18e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.25 on 6 degrees of freedom
Number of iterations to convergence: 11
Achieved convergence tolerance: 7.025e-06
## without derivative
nls(rate ~ Vm * conc/(K + conc),
data = L.minor, start = list(K = 1, Vm = 1),
trace = TRUE)
42634.85 (7.03e+00): par = (1 1)
31698.57 (9.47e+00): par = (389.6406 100.9883)
21357.27 (1.23e+01): par = (31.60848 46.75041)
18634.07 (1.03e+01): par = (27.26996 51.07311)
13989.24 (8.02e+00): par = (22.18639 59.6233)
7637.536 (5.64e+00): par = (18.35113 75.51722)
2003.136 (2.75e+00): par = (17.01927 100.4895)
234.3615 (6.31e-03): par = (17.1086 126.0567)
234.3533 (9.56e-04): par = (17.07326 126.021)
234.3531 (1.86e-04): par = (17.08017 126.0352)
234.3531 (3.61e-05): par = (17.07882 126.0324)
234.3531 (7.03e-06): par = (17.07909 126.0329)
Nonlinear regression model
model: rate ~ Vm * conc/(K + conc)
data: L.minor
K Vm
17.08 126.03
residual sum-of-squares: 234.4
Number of iterations to convergence: 11
Achieved convergence tolerance: 7.029e-06
MMfct2 <- deriv(~ Vm * conc / (K + conc),
c("K", "Vm"),
function(conc, K, Vm) {})
MMfct2
function (conc, K, Vm)
{
.expr1 <- Vm * conc
.expr2 <- K + conc
.value <- .expr1/.expr2
.grad <- array(0, c(length(.value), 2L), list(NULL, c("K",
"Vm")))
.grad[, "K"] <- -(.expr1/.expr2^2)
.grad[, "Vm"] <- conc/.expr2
attr(.value, "gradient") <- .grad
.value
}
MMfct2(1, 2, 3)
[1] 1
attr(,"gradient")
K Vm
[1,] -0.3333333 0.3333333
L.minor.mgr2 <- nls(
rate ~ MMfct2(conc, + K, Vm),
data = L.minor,
start = list(K = 20, Vm = 120))
summary(L.minor.mgr2)
Formula: rate ~ MMfct2(conc, +K, Vm)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 17.079 2.953 5.784 0.00117 **
Vm 126.033 7.173 17.570 2.18e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.25 on 6 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 8.148e-06
Ref: Venables와 Ripley(2002a, pp. 214–216)
L.minor.m3 <- nls(rate ~ conc/(K + conc),
data = L.minor,
algorithm = "plinear",
start = list(K = 20))
summary(L.minor.m3)
Formula: rate ~ conc/(K + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 17.079 2.953 5.784 0.00117 **
.lin 126.033 7.173 17.570 2.18e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.25 on 6 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 2.143e-06
## Fig. 4.1.
data(segreg)
plot(C ~ Temp, data = segreg,
xlab = Mean ~ temperature ~ (degree ~ F),
ylab = "Energy consumption (kWh)")
profRSS1 <- function(gamma) {
deviance(lm(C ~ pmax(0, Temp - gamma), data = segreg))
}
profRSS2 <- Vectorize(profRSS1, "gamma")
l <- lm( C ~ pmax(0, Temp - 40), data = segreg)
summary(l)
Call:
lm(formula = C ~ pmax(0, Temp - 40), data = segreg)
Residuals:
Min 1Q Median 3Q Max
-17.472 -2.155 -0.148 2.678 11.775
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 74.43931 1.16428 63.936 < 2e-16 ***
pmax(0, Temp - 40) 0.53635 0.06226 8.614 2.27e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.312 on 37 degrees of freedom
Multiple R-squared: 0.6673, Adjusted R-squared: 0.6583
F-statistic: 74.21 on 1 and 37 DF, p-value: 2.268e-10
plot(segreg$Temp, segreg$C)
lines(segreg$Temp, predict(l), col = 2)
## Fig. 4.2.
plot(profRSS2(Temp) ~ Temp,
data = segreg,
type = "l",
xlab = expression(gamma),
ylab = "Profile RSS")
## 좀더 부드럽게 만듬
plot(profRSS2(Temp) ~ Temp,
data = data.frame(Temp = seq(min(segreg$Temp),
max(segreg$Temp),
length.out = 1001)),
type = "l",
xlab = expression(gamma),
ylab = "Profile RSS")
RScompetition.m1 <- nls(
biomass ~ a/(1 + b * (x + c * z)),
data = RScompetition,
start = list(a = 20, b = 1, c = 1)
)
summary(RScompetition.m1)
Formula: biomass ~ a/(1 + b * (x + c * z))
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 25.9144 2.3393 11.078 1.42e-14 ***
b 0.1776 0.0355 5.004 8.68e-06 ***
c 1.1349 0.2964 3.829 0.000387 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.673 on 46 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 1.437e-06
virDensity <- with(RScompetition, x + coef(RScompetition.m1)[3] * z)
virDenVal <- seq(0, max(virDensity), length.out = 100)
biomassVal <- predict(RScompetition.m1, data.frame(x = virDenVal, z = 0))
## Fig. 4.3.
plot(biomassVal ~ virDenVal, type = "l",
ylab = "Biomass of sensitive biotype (g/plant)",
xlab = Virtual ~ density ~ (plants/m^2))
with(RScompetition, points(biomass ~ virDensity))
## Fig. 4.4.
par(pty = "s")
plot(Q ~ I, data = IQsig)
theta <- 0:360 * (pi/180)
lines(cos(theta), sin(theta))
IQsig.m1 <- nls(
~ ((I - I0)^2 - 2 * gamma * sin(phi) * (I - I0) * (Q - Q0) +
gamma * gamma * (Q - Q0)^2) - (rho * gamma * cos(phi))^2,
data = IQsig,
start = list(I0 = -0.005, gamma = 1, phi = -0.005, Q0 = -0.005, rho = 1))
summary(IQsig.m1)
Formula: 0 ~ ((I - I0)^2 - 2 * gamma * sin(phi) * (I - I0) * (Q - Q0) +
gamma * gamma * (Q - Q0)^2) - (rho * gamma * cos(phi))^2
Parameters:
Estimate Std. Error t value Pr(>|t|)
I0 -0.002259 0.002058 -1.097 0.278
gamma 1.007053 0.003719 270.764 < 2e-16 ***
phi -0.053540 0.004977 -10.758 5.02e-14 ***
Q0 -0.002481 0.002278 -1.089 0.282
rho 0.974295 0.002540 383.553 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01917 on 45 degrees of freedom
Number of iterations to convergence: 3
Achieved convergence tolerance: 4.699e-07
## Fig. 5.1.
plot(p ~ T, data = vapCO, log = "y",
xlab = "Temperature (K)",
ylab = "Pressure (Pa)")
vapCO.m1 <- nls(
log(p) ~ A - B/(C + T),
data = vapCO,
start = list(A = 10, B = 100, C = -10))
summary(vapCO.m1)
Formula: log(p) ~ A - B/(C + T)
Parameters:
Estimate Std. Error t value Pr(>|t|)
A 19.0923 0.2312 82.56 < 2e-16 ***
B 497.8246 27.7002 17.97 4.84e-10 ***
C -16.1516 1.5437 -10.46 2.19e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.09951 on 12 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 3.585e-06
## Fig. 5.2.
plot(p ~ T, data = vapCO, log = "y",
xlab = "Temperature (K)",
ylab = "Pressure (Pa)")
lines(vapCO$T, exp(fitted(vapCO.m1)), col = 2)
## Fig. 5.3.
plot(weight ~ conc, data = lettuce,
xlab = "Concentration (mg/l)",
ylab = "Biomass (g)", log = "x")
Warning in xy.coords(x, y, xlabel, ylabel, log) :
2 x values <= 0 omitted from logarithmic plot
## Fig. 5.4.
plot(fitted(vapCO.m1), residuals(vapCO.m1),
xlab = "Fitted Values", ylab = "Residuals")
abline(a = 0, b = 0)
## Fig. 5.5.
plot(rootl ~ conc, data = ryegrass,
xlab = "Concentration (mM)",
ylab = "Root length (cm)")
ryegrass.m1 <- lm(rootl ~ as.factor(conc), data = ryegrass)
summary(ryegrass.m1)
Call:
lm(formula = rootl ~ as.factor(conc), data = ryegrass)
Residuals:
Min 1Q Median 3Q Max
-1.08968 -0.17121 0.01577 0.16847 1.21865
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.74935 0.22535 34.388 < 2e-16 ***
as.factor(conc)0.94 -0.07606 0.39032 -0.195 0.84780
as.factor(conc)1.88 -1.33479 0.39032 -3.420 0.00327 **
as.factor(conc)3.75 -4.73466 0.39032 -12.130 8.53e-10 ***
as.factor(conc)7.5 -6.71542 0.39032 -17.205 3.45e-12 ***
as.factor(conc)15 -7.07018 0.39032 -18.114 1.50e-12 ***
as.factor(conc)30 -7.44601 0.39032 -19.077 6.47e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.552 on 17 degrees of freedom
Multiple R-squared: 0.9791, Adjusted R-squared: 0.9718
F-statistic: 133 on 6 and 17 DF, p-value: 2.488e-13
ryegrass.m2 <- nls(
rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e)))),
data = ryegrass,
start = list(b = 1, c = 0.6, d = 8, e = 3)
)
summary(ryegrass.m2)
Formula: rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e))))
Parameters:
Estimate Std. Error t value Pr(>|t|)
b 2.9822 0.4584 6.505 2.43e-06 ***
c 0.4814 0.2095 2.298 0.0325 *
d 7.7930 0.1897 41.073 < 2e-16 ***
e 3.0580 0.1858 16.456 4.31e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5196 on 20 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 6.33e-07
anova(ryegrass.m2, ryegrass.m1)
Analysis of Variance Table
Model 1: rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e))))
Model 2: rootl ~ as.factor(conc)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 20 5.4002
2 17 5.1799 3 0.22035 0.2411 0.8665
Q <- -2 * (logLik(ryegrass.m2) - logLik(ryegrass.m1))
df.Q <- df.residual(ryegrass.m2) - df.residual(ryegrass.m1)
1 - pchisq(Q, df.Q)
'log Lik.' 0.8012878 (df=5)
with(ryegrass, leveneTest(rootl, as.factor(conc)))
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 6 1.9266 0.1344
17
## Fig. 5.6.
plot(fitted(vapCO.m1), abs(residuals(vapCO.m1)),
xlab = "Fitted values", ylab = "Absolute residuals")
## Fig. 5.7.
plot(vapCO.m1)
## Fig. 5.8.
standardRes <- residuals(ryegrass.m2)/summary(ryegrass.m2)$sigma
par(pty = "s")
qqnorm(standardRes, main = "")
abline(a = 0, b = 1)
shapiro.test(standardRes)
Shapiro-Wilk normality test
data: standardRes
W = 0.98232, p-value = 0.9345
## Fig. 5.9.
plot(residuals(vapCO.m1), c(residuals(vapCO.m1)[-1], NA),
xlab = "Residuals", ylab = "Lagged residuals")
## Fig. 6.1.
plot(RGR ~ Day, data = RGRcurve,
xlab = "Time (days)", ylab = "Relative growth rate (%)")
y0_start = mean(RGRcurve[RGRcurve$Day == 0, "RGR"])
y1 = mean(RGRcurve[RGRcurve$Day == 1, "RGR"])
b_start = 1 / log(y1/y0_start)
RGRcurve.m2 <- gnls(
RGR ~ y0 * exp(Day/b),
data = RGRcurve,
start = list(y0 = y0_start, b = b_start),
weights = varPower()
)
summary(RGRcurve.m2)
Generalized nonlinear least squares fit
Model: RGR ~ y0 * exp(Day/b)
Data: RGRcurve
Variance function:
Structure: Power of variance covariate
Formula: ~fitted(.)
Parameter estimates:
power
1.012821
Coefficients:
Correlation:
y0
b 0.797
Standardized residuals:
Min Q1 Med Q3 Max
-1.8322566 -0.6356287 -0.1612768 0.6476686 1.8176689
Residual standard error: 0.1672069
Degrees of freedom: 41 total; 39 residual
RGRcurve.m2g <- gnls(
RGR ~ y0 * exp(Day/b),
data = RGRcurve,
start = list(y0 = y0_start, b = b_start)
# weights = varPower() # without weights
)
summary(RGRcurve.m2g )
Generalized nonlinear least squares fit
Model: RGR ~ y0 * exp(Day/b)
Data: RGRcurve
Coefficients:
Correlation:
y0
b 0.962
Standardized residuals:
Min Q1 Med Q3 Max
-2.4946822 -0.5595546 -0.1065972 0.3415913 2.3870757
Residual standard error: 0.1009825
Degrees of freedom: 41 total; 39 residual
AIC(RGRcurve.m2, RGRcurve.m2g )
## Fig. 6.2.
plot(recruits ~ spawners,
data = sockeye[-12,],
xlab = "Number of spawners (thousands)",
ylab = "Number of recruits (thousands)")
일단 nls를 이용하여 fitting을 한 후
sockeye.m1 <- nls(
recruits ~ beta1 * spawners * exp(-beta2 * spawners),
data = sockeye[-12, ],
start = list(beta1 = 2, beta2 = 0.001)
)
summary(sockeye.m1)
Formula: recruits ~ beta1 * spawners * exp(-beta2 * spawners)
Parameters:
Estimate Std. Error t value Pr(>|t|)
beta1 3.9495589 1.0154151 3.890 0.000658 ***
beta2 0.0008526 0.0003570 2.388 0.024812 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 566.1 on 25 degrees of freedom
Number of iterations to convergence: 4
Achieved convergence tolerance: 3.796e-07
## Fig. 6.3.
plot(fitted(sockeye.m1), abs(residuals(sockeye.m1)),
xlab = "Fitted values", ylab = "Absolute residuals")
boxcox.nls를 이용하여 양쪽을 변화시키며 피팅을 한다.
## Fig. 6.4.
sockeye.m2 <- boxcox.nls(sockeye.m1)
bcSummary(sockeye.m2)
Estimated lambda: -0.2
Confidence interval for lambda: [-0.89, 0.52]
coef(summary(sockeye.m1))
Estimate Std. Error t value Pr(>|t|)
beta1 3.9495588725 1.0154150955 3.889600 0.000657508
beta2 0.0008525523 0.0003570088 2.388043 0.024811780
coef(summary(sockeye.m2))
Estimate Std. Error t value Pr(>|t|)
beta1 3.7767882170 0.6941022768 5.441256 1.195243e-05
beta2 0.0009539343 0.0003061565 3.115839 4.563494e-03
이 내용은 정확히 이해하지는 못했음.
vcov(sockeye.m1)
beta1 beta2
beta1 1.0310678161 3.436055e-04
beta2 0.0003436055 1.274553e-07
sandwich(sockeye.m1)
beta1 beta2
beta1 0.6082739015 2.458975e-04
beta2 0.0002458975 1.162302e-07
coeftest(sockeye.m1)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
beta1 3.94955887 1.01541510 3.8896 0.0006575 ***
beta2 0.00085255 0.00035701 2.3880 0.0248118 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
coeftest(sockeye.m1, vcov = sandwich)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
beta1 3.94955887 0.77991916 5.0641 3.158e-05 ***
beta2 0.00085255 0.00034093 2.5007 0.01931 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp1
## Fig. 6.5.
plot(Nremaining ~ time, data = exp1,
xlab = "Time (years)", ylab = "Nitroge content (%)")
exp1.m1 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp1)
exp1.m2 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp1,
weights = norep/(stdev * stdev))
weights(exp1.m2)
[1] 0.4403928 0.4403928 0.1096776 0.1430179 1.1291355 0.0867694 1.8896447 0.3667661
coef(summary(exp1.m1))
Estimate Std. Error t value Pr(>|t|)
a1 49.6945735 6.1714388 8.0523481 0.001291449
a2 0.2182533 0.3162611 0.6901049 0.528080839
b1 32.2059455 6.1631160 5.2255946 0.006402929
b2 -3.1242900 0.4881789 -6.3998879 0.003061147
coef(summary(exp1.m2))
Estimate Std. Error t value Pr(>|t|)
a1 48.95580529 5.8454395 8.37504275 0.001111764
a2 -0.02288415 0.2953566 -0.07747973 0.941962763
b1 30.67388498 6.4550371 4.75193009 0.008958175
b2 -3.15057479 0.5205589 -6.05229290 0.003760909
exp2
exp2.m1 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp2)
exp2.m2 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp2,
weights = norep/(stdev^2))
## Fig. 6.6.
plot(Nremaining ~ time, data = exp2,
xlab = "Time (years)", ylab = "Nitrogen content (%)")
timeVal <- with(exp2, seq(min(time), max(time), length.out = 101))
lines(timeVal, predict(exp2.m1, newdata = data.frame(time = timeVal)),
lty = 2, col = 2)
lines(timeVal, predict(exp2.m2, newdata = data.frame(time = timeVal)),
lty = 3, col = 3)
coef(summary(exp2.m1))
Estimate Std. Error t value Pr(>|t|)
a1 37.1828528 7.3946625 5.0283367 0.007342087
a2 0.3984263 0.4676840 0.8519136 0.442262753
b1 56.9940958 7.2995778 7.8078620 0.001451984
b2 -2.4593036 0.2429193 -10.1239548 0.000535818
coef(summary(exp2.m2))
Estimate Std. Error t value Pr(>|t|)
a1 28.9158392 8.3745109 3.4528392 2.598700e-02
a2 0.3128281 0.8122243 0.3851499 7.197304e-01
b1 59.7777230 3.7595533 15.9002197 9.144765e-05
b2 -2.3567707 0.1311759 -17.9664898 5.641332e-05
summary(exp1.m1)
Formula: Nremaining ~ SSbiexp(time, a1, a2, b1, b2)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a1 49.6946 6.1714 8.052 0.00129 **
a2 0.2183 0.3163 0.690 0.52808
b1 32.2059 6.1631 5.226 0.00640 **
b2 -3.1243 0.4882 -6.400 0.00306 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.844 on 4 degrees of freedom
Number of iterations to convergence: 0
Achieved convergence tolerance: 5.364e-06
summary(exp1.m2)
Formula: Nremaining ~ SSbiexp(time, a1, a2, b1, b2)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a1 48.95581 5.84544 8.375 0.00111 **
a2 -0.02288 0.29536 -0.077 0.94196
b1 30.67388 6.45504 4.752 0.00896 **
b2 -3.15057 0.52056 -6.052 0.00376 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.785 on 4 degrees of freedom
Number of iterations to convergence: 9
Achieved convergence tolerance: 9.037e-06
L.minor.m1 <- update(L.minor.m1, trace = FALSE)
L.minor.m1pro <- profile(L.minor.m1)
## Fig. 7.1.
plot(L.minor.m1pro)
## Fig. 7.2.
plot(L.minor.m1pro, absVal = FALSE)
confint(L.minor.m1)
Waiting for profiling to be done...
2.5% 97.5%
K 11.45142 24.86684
Vm 110.64232 143.75648
set.seed(222)
L.minor.m1boot <- nlsBoot(L.minor.m1)
summary(L.minor.m1boot)
------
Bootstrap statistics
Estimate Std. error
K 17.37534 2.730157
Vm 126.64954 6.522141
------
Median of bootstrap estimates and percentile confidence intervals
Median 2.5% 97.5%
K 17.25307 12.30746 23.20204
Vm 126.08120 115.00875 139.46932
## Fig. 7.3.
par(mfrow = c(1, 2))
qqnorm(L.minor.m1boot$coefboot[, 1], main = "K")
qqnorm(L.minor.m1boot$coefboot[, 2], main = "Vm")
confint2(L.minor.m1)
2.5 % 97.5 %
K 9.853828 24.30416
Vm 108.480813 143.58470
confint2(L.minor.m1, level = 0.99)
0.5 % 99.5 %
K 6.131815 28.02617
Vm 99.439005 152.62651
deltaMethod(L.minor.m1, "Vm/(4*K)")
Estimate SE 2.5 % 97.5 %
Vm/(4 * K) 1.8449 0.2363 1.3817 2.308
secalonic
## Fig. 7.4.
plot(rootl ~ dose, data = secalonic,
xlab = "Dose (mM)", ylab = "Root length (cm)")
secalonic.m1 <- nls(rootl ~ SSfpl(dose, a, b, c, d), data = secalonic)
summary(secalonic.m1)
Formula: rootl ~ SSfpl(dose, a, b, c, d)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 6.053612 0.395467 15.308 0.000606 ***
b 0.353944 0.194089 1.824 0.165722
c 0.075188 0.005911 12.721 0.001048 **
d 0.029350 0.006621 4.433 0.021333 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2028 on 3 degrees of freedom
Number of iterations to convergence: 0
Achieved convergence tolerance: 8.102e-06
secalonic.m2 <- nls(rootl ~ SSlogis(dose, a, c, d), data = secalonic)
summary(secalonic.m2)
Formula: rootl ~ SSlogis(dose, a, c, d)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 6.344853 0.564690 11.236 0.000357 ***
c 0.076959 0.009397 8.190 0.001211 **
d -0.035740 0.007244 -4.933 0.007854 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2546 on 4 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 4.999e-07
anova(secalonic.m2, secalonic.m1)
Analysis of Variance Table
Model 1: rootl ~ SSlogis(dose, a, c, d)
Model 2: rootl ~ SSfpl(dose, a, b, c, d)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 4 0.25924
2 3 0.12341 1 0.13582 3.3016 0.1668
M.merluccius
plot(num.fish ~ spawn.biomass, data = M.merluccius)
4가지 조금씩 다른 모델을 비교
M.merluccius.bh <- nls(
num.fish ~ spawn.biomass * alpha/(1 + spawn.biomass/k),
data = M.merluccius, start = list(alpha = 5, k = 50)
)
M.merluccius.de <- nls(
num.fish ~ spawn.biomass * alpha * (1 - c * spawn.biomass/k)^(1/c),
data = M.merluccius, start = list(alpha = 4.4, k = 106, c = 0.86)
)
M.merluccius.ri <- nls(
num.fish ~ spawn.biomass * alpha * exp(-spawn.biomass/k),
data = M.merluccius, start = list(alpha = 5, k = 50)
)
M.merluccius.sh <- nls(
num.fish ~ spawn.biomass * alpha/(1 + (spawn.biomass/k)^c),
data = M.merluccius, start = list(alpha = 3.87, k = 61.72, c = 2.25),
control = nls.control(maxiter = 100)
)
Residual standard error를 계산하면,
summary(M.merluccius.bh)$sigma
[1] 14.69956
summary(M.merluccius.de)$sigma
[1] 14.78346
summary(M.merluccius.ri)$sigma
[1] 14.30425
summary(M.merluccius.sh)$sigma
[1] 14.56758
AIC를 비교해보면,
AIC(M.merluccius.bh)
[1] 127.0562
AIC(M.merluccius.de)
[1] 128.0263
AIC(M.merluccius.ri)
[1] 126.2383
AIC(M.merluccius.sh)
[1] 127.585
Burnham and Anderson (2002, pp. 70–72)은 한 모델이 다른 모델에 비하여 절대적으로 우위에 있기 위해서는 AIC의 차이가 10 이상의 차이가 필요하다는 경험적 법칙을 제시하였다. 따라서 4 모델 사이에서 큰 차이는 없다.
Puromycin
## Fig. 8.1.
xyplot(rate ~ conc | state, data = Puromycin,
xlab = "Substrate concentration (ppm)",
ylab = "Reaction rates\n(counts/min/min)")
Puromycin.m1 <- nls(
rate ~ Vm[state] * conc/(K[state] + conc),
data = Puromycin,
start = list(K = c(0.1, 0.1), Vm = c(200, 200))
)
summary(Puromycin.m1)
Formula: rate ~ Vm[state] * conc/(K[state] + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K1 6.412e-02 7.877e-03 8.141 1.29e-07 ***
K2 4.771e-02 8.281e-03 5.761 1.50e-05 ***
Vm1 2.127e+02 6.608e+00 32.185 < 2e-16 ***
Vm2 1.603e+02 6.896e+00 23.242 2.04e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.4 on 19 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 1.568e-06
Puromycin.m2 <- gnls(
rate ~ Vm * conc/(K + conc),
data = Puromycin,
start = list(Vm = c(200, 200), K = c(0.1, 0.1)),
params = list(Vm ~ state - 1, K ~ state - 1)
)
summary(Puromycin.m2)
Generalized nonlinear least squares fit
Model: rate ~ Vm * conc/(K + conc)
Data: Puromycin
Coefficients:
Correlation:
Vm.sttt Vm.sttn K.sttt
Vm.stateuntreated 0.000
K.statetreated 0.765 0.000
K.stateuntreated 0.000 0.777 0.000
Standardized residuals:
Min Q1 Med Q3 Max
-1.32632770 -0.52742178 -0.06890298 0.51295804 2.44557119
Residual standard error: 10.40003
Degrees of freedom: 23 total; 19 residual
Puromycin.m3 <- nlsList(
rate ~ SSmicmen(conc, Vm, K) | state,
data = Puromycin
)
summary(Puromycin.m3)
Call:
Model: rate ~ SSmicmen(conc, Vm, K) | state
Data: Puromycin
Coefficients:
Vm
Estimate Std. Error t value Pr(>|t|)
treated 212.6837 6.608092 32.18534 3.241160e-11
untreated 160.2801 6.896019 23.24241 1.384618e-09
K
Estimate Std. Error t value Pr(>|t|)
treated 0.06412123 0.007876786 8.140532 1.565136e-05
untreated 0.04770829 0.008281171 5.761056 1.727047e-04
Residual standard error: 10.40003 on 19 degrees of freedom
Puromycin2 <- groupedData(rate ~ conc | state, data = Puromycin)
Puromycin.m4 <- nlsList(
rate ~ SSmicmen(conc, a, b), data = Puromycin2
)
summary(Puromycin.m4)
Call:
Model: rate ~ SSmicmen(conc, a, b) | state
Data: Puromycin2
Coefficients:
a
Estimate Std. Error t value Pr(>|t|)
untreated 160.2801 6.896019 23.24241 1.384618e-09
treated 212.6837 6.608092 32.18534 3.241160e-11
b
Estimate Std. Error t value Pr(>|t|)
untreated 0.04770829 0.008281171 5.761056 1.727047e-04
treated 0.06412123 0.007876786 8.140532 1.565136e-05
Residual standard error: 10.40003 on 19 degrees of freedom
Puromycin.m5 <- nls(
rate ~ Vm * conc/(K + conc), data = Puromycin,
start = list(K = 0.1, Vm = 200))
summary(Puromycin.m5)
Formula: rate ~ Vm * conc/(K + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 0.06039 0.01077 5.608 1.45e-05 ***
Vm 190.80620 8.76458 21.770 6.84e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 18.61 on 21 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 6.259e-06
anova(Puromycin.m5, Puromycin.m1)
Analysis of Variance Table
Model 1: rate ~ Vm * conc/(K + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 21 7276.5
2 19 2055.1 2 5221.5 24.138 6.075e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AIC(Puromycin.m5, Puromycin.m1)
Puromycin.m6 <- nls(
rate ~ Vm * conc/(K[state] + conc),
data = Puromycin,
start = list(K = c(0.1, 0.1), Vm = 200))
summary(Puromycin.m6)
Formula: rate ~ Vm * conc/(K[state] + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K1 0.04801 0.00874 5.493 2.24e-05 ***
K2 0.08865 0.01448 6.122 5.54e-06 ***
Vm 194.10328 7.40464 26.214 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.52 on 20 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 8.957e-06
anova(Puromycin.m6, Puromycin.m1)
Analysis of Variance Table
Model 1: rate ~ Vm * conc/(K[state] + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 20 4815.4
2 19 2055.1 1 2760.3 25.521 7.082e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Puromycin.m7 <- nls(
rate ~ Vm[state] * conc/(K + conc),
data = Puromycin, start = list(K = 0.1, Vm = c(200, 200))
)
summary(Puromycin.m7)
Formula: rate ~ Vm[state] * conc/(K + conc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
K 0.05797 0.00591 9.809 4.37e-09 ***
Vm1 208.63004 5.80399 35.946 < 2e-16 ***
Vm2 166.60408 5.80743 28.688 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.59 on 20 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 1.364e-06
anova(Puromycin.m7, Puromycin.m1)
Analysis of Variance Table
Model 1: rate ~ Vm[state] * conc/(K + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 20 2240.9
2 19 2055.1 1 185.84 1.7182 0.2056
AIC로 비교해 보아도 Puromycin.m7이 더욱 좋다.
AIC(Puromycin.m1, Puromycin.m5, Puromycin.m6, Puromycin.m7)
G.aparine
## Fig. 8.2.
xyplot(drymatter ~ dose | as.factor(treatment),
data = G.aparine,
xlab = "Dose (g/ha)", ylab = "Dry matter (mg/pot)")
이미 데이터에서 treatment 0가 없어졌음. 아마도 데이터를 업데이트 한 것으로 추측됨.
G.aparine$treatment2 <- factor(G.aparine$treatment)
#levels(G.aparine$treatment2) <- c("0", "1", "2")
G.aparine
G.aparine.m1 <- nls(
drymatter ~ c[treatment2] + (d - c[treatment2])/(1 + exp(b[treatment2] * (log(dose) - log(e[treatment2])))),
data = G.aparine,
start = list(b = c(2, 2), c = c(500, 100), d = 1000, e = c(50, 100)))
summary(G.aparine.m1)
Formula: drymatter ~ c[treatment2] + (d - c[treatment2])/(1 + exp(b[treatment2] *
(log(dose) - log(e[treatment2]))))
Parameters:
Estimate Std. Error t value Pr(>|t|)
b1 1.6129 0.3462 4.659 5.33e-06 ***
b2 1.7510 0.2303 7.603 7.07e-13 ***
c1 509.5035 23.4379 21.738 < 2e-16 ***
c2 151.9164 27.2114 5.583 6.56e-08 ***
d 984.8883 12.7957 76.970 < 2e-16 ***
e1 50.8000 7.8465 6.474 5.58e-10 ***
e2 93.4466 8.1626 11.448 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 111.4 on 233 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 9.945e-06
concValues <- with(Puromycin, seq(min(conc), max(conc, length.out = 11)))
concValues
[1] 0.02 1.02 2.02 3.02 4.02 5.02 6.02 7.02 8.02 9.02 10.02
stateVal1 <- levels(Puromycin$state)
stateVal1
[1] "treated" "untreated"
csValues1 <- expand.grid(conc = concValues, state = stateVal1)
csValues1
predict(Puromycin.m7, newdata = csValues1)
[1] 53.51423 197.41021 202.80963 204.70062 205.66419 206.24825 206.64012
[8] 206.92127 207.13280 207.29773 207.42993 42.73445 157.64435 161.95611
[15] 163.46619 164.23566 164.70207 165.01500 165.23951 165.40844 165.54015
[22] 165.64572
stateVal2 <- factor("untreated", levels = c("treated", "untreated"))
stateVal2
[1] untreated
Levels: treated untreated
csValues2 <- data.frame(conc = concValues, state = stateVal2)
csValues2
predict(Puromycin.m1, newdata = csValues2)
[1] 47.3444 153.1183 156.5819 157.7874 158.4002 158.7711 159.0198 159.1981
[9] 159.3322 159.4367 159.5205
vinclozolin
## Fig. 8.3.
xyplot(effect ~ conc | exper, data = vinclozolin,
xlab = expression(paste("Concentration (", mu, "M)")),
ylab = "Luminescence (LU")
LL3.formula <- effect ~ d/(1 + exp(b * (log(conc) - log(e))))
vinclozolin.e1.m <- nls(
LL3.formula, data = vinclozolin, subset = exper == 10509,
start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e2.m <- nls(
LL3.formula, data = vinclozolin, subset = exper == 10821,
start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e3.m <- nls(
LL3.formula, data = vinclozolin, subset = exper == 10828,
start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e4.m <- nls(
LL3.formula, data = vinclozolin, subset = exper == 10904,
start = list(b = 1, d = 2700, e = 0.03))
vinclozolin.e5.m <- nls(
LL3.formula, data = vinclozolin, subset = exper == 11023,
start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e6.m <- nls(
LL3.formula, data = vinclozolin, subset = exper == 11106,
start = list(b = 0.5, d = 2600, e = 0.02))
## Fig. 8.4.
#par(pty = "s")
plot(effect ~ conc, data = vinclozolin,
pch = as.numeric(exper),
log = "x",
xlim = c(1e-04, 10),
xlab = expression(paste("Concentration(", mu, "M)")),
ylab = "Luminescence(LU)")
concVec <- exp(seq(log(1e-04), log(10), length.out = 50))
lines(concVec, predict(vinclozolin.e1.m, data.frame(conc = concVec)), lty = 2)
lines(concVec, predict(vinclozolin.e2.m, data.frame(conc = concVec)), lty = 3)
lines(concVec, predict(vinclozolin.e3.m, data.frame(conc = concVec)), lty = 4)
lines(concVec, predict(vinclozolin.e4.m, data.frame(conc = concVec)), lty = 5)
lines(concVec, predict(vinclozolin.e5.m, data.frame(conc = concVec)), lty = 6)
lines(concVec, predict(vinclozolin.e6.m, data.frame(conc = concVec)), lty = "3313")
vinclozolin.m1 <- nlme(
effect ~ d/(1 + exp(b * (log(conc) - log(e)))),
data = vinclozolin,
fixed = list(b ~ 1, d ~ 1, e ~ 1),
random = d ~ 1 | exper,
start = c(1, 1000, 0.1))
summary(vinclozolin.m1)
Nonlinear mixed-effects model fit by maximum likelihood
Model: effect ~ d/(1 + exp(b * (log(conc) - log(e))))
Data: vinclozolin
Random effects:
Formula: d ~ 1 | exper
d Residual
StdDev: 695.3775 150.9623
Fixed effects: list(b ~ 1, d ~ 1, e ~ 1)
Correlation:
b d
d -0.051
e 0.529 -0.145
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.83842036 -0.58871949 -0.06161384 0.54652838 2.07902663
Number of Observations: 53
Number of Groups: 6
## Fig. 8.5.
#par(pty = "s")
plot(effect ~ conc, data = vinclozolin,
pch = as.numeric(exper),
log = "x",
xlim = c(1e-04, 10),
xlab = expression(paste("Concentration(", mu, "M)")),
ylab = "Luminescence(LU)")
concVec <- exp(seq(log(1e-04), log(10), length.out = 50))
lines(concVec, predict(vinclozolin.e1.m, data.frame(conc = concVec)), lty = 2)
lines(concVec, predict(vinclozolin.e2.m, data.frame(conc = concVec)), lty = 3)
lines(concVec, predict(vinclozolin.e3.m, data.frame(conc = concVec)), lty = 4)
lines(concVec, predict(vinclozolin.e4.m, data.frame(conc = concVec)), lty = 5)
lines(concVec, predict(vinclozolin.e5.m, data.frame(conc = concVec)), lty = 6)
lines(concVec, predict(vinclozolin.e6.m, data.frame(conc = concVec)), lty = "3313")
lines(concVec,
predict(vinclozolin.m1, newdata = data.frame(conc = concVec), level = 0),
lty = 1, lwd = 3)